Phase transitions of a single polyelectrolyte in a poor solvent with explicit counterions 
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Abstract 

Conformational properties of a single flexible polyelectrolyte chain in a poor solvent are studied using constant temperature 
molecular dynamics simulation. The effects of counterions are explicitly taken in to account. Structural properties of various 
phases and the transition between these phases are studied by tracking the values of asphericity, radius of gyration, fraction 
of condensed counterions, number of non-bonded neighbours and Coulomb interaction energies. From our simulations, we 
find strong evidence for a first-order phase transition from extended to collapsed phase consistent with earlier theoretical 
predictions. We also identify a continuous phase transition associated with the condensation of counterions and estimate the 
critical exponents associated with the transition. Finally, we argue that previous suggestions of existence of an independent 
intermediate phase between extended and collapsed phases is only a finite size effect. 

PACS numbers: 82.35.Rs,36.20.Ey,64.70.km 



I. INTRODUCTION 

Polyelectrolytes are polymers, which in a polar sol- 
vent, release counterions into the solution, making the 
polymer backbone charged. Examples of polyelectrolytes 
include sulfonated polystyrene, polymethacrylic acid, 
DNA, RNA and proteins. Though the system is overall 
charge neutral, long range Coulomb interactions make 
the behaviour of polyelectrolytes considerably different 
from that of neutral polymers ^ . The static and dy- 
namic properties of polyelectrolytes depend on the nature 
of the solvent, linear charge density of the polymer back 
bone, valency of the counterions, temperature, salt con- 
centration, hydrodynamic interactions and the bending 
rigidity 

A polyelectrolyte chain in a poor solvent undergoes 
a series of phase transitions as the linear charge den- 
sity of the chain is varied keeping temperature fixed. 
The rich phase diagram resulting from a competition 
between the long range Coulomb interaction and the 
short range excluded volume interaction has been stud- 
ied theoretically numerically [gI-IsI. [i^ - [15| and 
experimentally |16| - |22| . At very low charge densities, the 
polyelectrolyte chain is in a collapsed phase, with the 
counterions uniformly distributed throughout the solu- 
tion [ll] . On increasing the charge density, the polyelec- 
trolyte chain makes a transition into the pearl-necklace 
phase in which the chain has two (dumbbell) or more 
globules connected through a string of monomers [4|-[7[. 
This instability of the charged collapsed phase is similar 
to the Rayleigh instability of a charged droplet 23, 24]. 
While the pearl-necklace phase has been observed in nu- 
merical simulatio ns |T^ , [isl [25j , the experimental status 
is unclear [l^, [l^, |21|. Further increase in charge density 
decreases the number and size of the globules and the 



•Electronic address: lanoopiSimsc.res.in 



i.in| 



•f Electronic address: ,rrajesh@imsc.res.in| 



polyelectrolyte chain becomes extended. When the lin- 
ear charge density exceeds a threshold value, the electro- 
static energy dominates over the thermal fluctuations and 
the counterions condense on the polyelectrolyte chain . 
The condensed counterions form dipoles with the poly- 
mer monomers and the attraction among the dipoles 
leads to the collapse of the polyelectrolyte chain [12| . Re- 
cent simulations also show that, in extreme poor solvent 
conditions, the polyelectrolyte chain can make a direct 
transition from the initial collapsed phase to the final 
condensed collapsed phase with out encountering some 
or even all the intermediate phases [Tsj . 

Certain features of the phase diagram remain poorly 
understood. Recent simulations argue for the existence 
of an intermediate phase between the extended and the 
condensed collapsed phases . Referred to as the 

sausage phase, this phase is defined as a collapsed phase 
in which the shape of the collapsed polymer becomes as- 
pherical, having a non-zero mean asphericity. It is not 
clear whether this phase is just a finite size effect. 

In addition, the nature of the counterion condensation 
is not well understood. By studying the condensation of 
counterions on a cylinder in three dimensions, and on a 
disc in two dimensions, it was shown that condensation 
is a second order phase transition [26'-'28'|. However, to 
the best of our knowledge, the corresponding question 
has not been addressed for a three dimensional polyelec- 
trolyte chain system. 

In this paper, we simulate polyelectrolyte chains of dif- 
ferent lengths and varying charge densities. We argue 
that the sausage phase does not exist and is an arte- 
fact of studying very small chains. Our simulations also 
show strong evidence that the counterion condensation 
is a second order transition accompanied by a divergence 
in the fluctuations of the number of non-bonded nearest 
neighbours of a monomer. In addition, we study aspects 
of the extended to collapsed transition. 
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II. MODEL AND SIMULATION METHOD 

We model the polyelectrolyte chain as spherical 
beads (monomers) connected through springs where each 
monomer carries a charge qe. Counterions are monova- 
lent and are modelled as spherical beads, each carrying a 
charge —qe. The polyelectrolyte chain and the counteri- 
ons are assumed to be in a medium of uniform dielectric 
constant e. The potential energy due to the pair of par- 
ticles i and j consists of three interactions. 

Coulomb interaction: The electrostatic energy is given 

by 
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Zq'e 



(1) 



where Z = —1 for monomer-counterion pairs and Z = 1 
otherwise, and Vij is the distance between particle i and 
j- 

Excluded volume interaction: The excluded volume in- 
teractions are modelled by the Lennard- Jones potential, 
which for two particles at a distance rij , is given by 
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where is the minimum of the potential and a is the in- 
ter particle distance at which the potential becomes zero. 
We use reduced units, in which the energy and length 
scales are specified in units of (counterion-counterion) 
and a respectively. The depth of the attractive po- 
tential is chosen as 1.0 for monomer-counterion and 
counterion-counterion pairs and 2.0 for monomer pairs, 
while a is set to 1.0 for all pairs. We use the shifted 
Lennard- Jones potential in which ULj{rij) is set to zero 
beyond a cut off distance r^. The value of Tc equals 
1.0 for monomer-counterion and counterion-counterion 
pairs and 2.5 for monomer-monomer pairs. With this 
choice of the parameters, the excluded volume interac- 
tion is purely repulsive for all the pairs other than the 
monomer-monomer pairs. The effective short range at- 
traction among the monomers mimics poor solvent con- 
ditions. Other ways of realising poor solvent conditions 
may be found in Ref. f25| . 

Bond stretching interaction: The bond stretching en- 
ergy for pairs in the polymer that are connected directly 
through springs is given by 



(3) 



where k is the spring constant and b is the equilibrium 
bond length. The values of k and b are taken as 500 and 
1.12 respectively. This value of b is close to the minimum 
of Lennard- Jones potential. 

The relative strength of the electrostatic interaction is 
parametrised by a dimensionless quantity A: 
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where is the Bjerrum length |29|, the length scale 
below which electrostatic interaction dominates thermal 
fluctuations. 
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where fc^ is the Boltzmann constant and T is tempera- 
ture. The thermal fluctuations dominate over the elec- 
trostatic interactions for very low values of A and the 
counterions will be uniformly distributed in the solution. 
When A is of order one, the electrostatic interaction en- 
ergy is comparable to the thermal energy and the counte- 
rions begin to undergo Manning condensation 9| . In our 
simulation, we vary A from 0.055 to 14.29. 

The equations of motion are integrated in time 
using the molecular dynamics simulation package 
LAMMPS [H The simulations are carried out at 
constant temperature (T=1.0), maintained through a 
Nose-Hoover thermostat (coupling constant — 0.1) [H, 
[33} . The system is placed in a cubic box with periodic 
boundary conditions. At the start of the simulations, 
the conflguration of the chain is randomly chosen and 
the monovalent counterions are distributed uniformly 
throughout the volume such that the charge neutrality is 
achieved. In our simulations, the length N of the chain is 
varied from 50 to 400, keeping the overall particle density 
of the system fixed at a constant value, 7.23 x 10^^ par- 
ticles per . At this density, there is no direct contact 
between the polyelectrolyte chain and any of its peri- 
odic images. We use the particle-particle/particle- mesh 
(PPPM) technique Q to evaluate the energy and forces 
due to the long range Coulomb interactions. The time 
step of the integration is chosen as 0.001. Our simulation 
runs are divided into an equilibration run (5 X 10^ steps), 
followed by a production run (6x10^ steps). All the data 
shown in the plots are measurements over only the pro- 
duction run (1.2 x 10^ configurations). The errors in the 
plots are estimated by dividing the production into five 
statistically independent blocks and measuring the stan- 
dard deviation of the mean values of the observable in 
each block ,35.] . 



III. RESULTS AND DISCUSSION 

A. Configurations of the polyelectrolyte chain 

Snap shots of the polyelectrolyte configuration for dif- 
ferent values of A are shown in Fig. [TJ For very small 
values of A, the electrostatic interactions can be ignored, 
and the polymer exists in a collapsed phase [Fig.[lja)]. As 
A is increased, the globule breaks up into two [Fig. [Ifb)] 
or more smaller globules [Fig. [IJc)]. This pearl- necklace 
configuration becomes extended on further increasing A 
[Fig.[ljd)]. Counterion condensation is initiated and for 
sufficiently large A, the polymer undergoes a collapse 
transition to form a sausage phase [iBj [Fig. (He)] or a 
spherical globule [Fig. [Iff)]. 
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(e) A = 6.752 (f) A = 14.285 



FIG. 1: Snapshots of the polyelectrolyte chain configurations 
for different values of A. For the sake of clarity, counterions 
are not shown and figures are not in the same scale, (a) Glob- 
ular phase, (b) Dumbbell phase, (c) Pearl-necklace phase, (d) 
Extended phase, (e) Sausage phase and (f) Globular phase. 
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B. Distribution of counterions 



We consider a counterion to be condensed if its dis- 
tance from any monomer is less than 2£b 12]. Let Nc be 
the number of condensed counterions and N the number 
of monomers. The mean fraction of condensed counte- 
rions (Nc/N) as a function of A is shown in Fig. [5] for 
different N. With the above definition of a condensed 
couterion, it is observed that the counterion condensation 
starts slightly below A = 1, consistent with the findings 
in Ref . [3^ . Note that for the case of a infinitely long and 
uniformly charged cylinder, the classic Manning conden- 
sation occurs at ^ = 1. 

For a given value of A, the fraction of condensed coun- 
terions increases with increasing N and reaches a limiting 
value. A similar result was obtained for the good solvent 
case as well JJj . The dependence of the fraction of con- 
densed counterions on N is closely linked to the effective 
size or morphology of the polyelectrolyte chain. We show 
that the longer chains have smaller relative size for values 
of A at which the counterion condensation occurs. The 
dependence of the relative size of the chain on N can be 
studied by calculating the radius of gyration Rg, defined 



where fi is the position of the i particle and rem is 
the centre of mass of the chain. For fixed A, Rg/N, a 
measure of the relative extension, decreases with N in 
the region A ^ 0.89 where the condensation occurs (see 
Fig.©. 

These two observations, the effective size being smaller 
for larger chains, and longer chains having a larger frac- 
tion of condensed counterions (see Fig. [5]), are related. 
We argue that a polymer with a smaller effective size has 
more condensed counterions. In the inset of Fig. [51 we 
show that increasing emm (depth of the Lennard- Jones 
potential of monomer-monomer pairs) , keeping other pa- 
rameters fixed, results in increased condensation. It is 
clear that increasing emm can only result in the effective 
size of the polyelectrolyte chain becoming smaller. The 
increased condensation may be due to counterions expe- 
riencing lower electrostatic potential for a more compact 
chain. 

We also observe that Rg/N shows a jump around 
Ac ~ 6.25. The jump is more pronounced for longer 
chains. For A > 6.25, radius of gyration Rg scales as 
N^/^, indicating that the chain is in a collapsed phase. 
We also observe a small kink in the fraction of condensed 
ions (see Fig.[21) around the same value of A where Rg/N 
shows a jump. 
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FIG. 3: Ratio of the radius of gyration Rg to the number of 
monomers A'^ as a function of A. Inset: The data for RgjN^^^ 
collapse for different A*' for A > 6.25. 



C. Number of non-bonded neighbours and the 
condensation transition 

The number of non-bonded nearest neighbours of a 
monomer is a useful order parameter in studying the 
collapse transition of a neutral polymer [33|. We study 
its behaviour for the polyelectrolyte chain. For a given 
monomer, a non-bonded neighbour is defined as any 
monomer/counterion that is not connected to it by a 
bond and within a distance b. The variation of the 
mean number of non bonded neighbours per monomer 
(rib) with A is shown in Fig. S) {rib) decreases when 
the polyelectrolyte chain goes from the initial collapsed 
phase to the extended phase. It has the minimum around 
A' « 0.89 where the chain is most extended. Close to 
Ac ~ 6.25, (rift) has a jump with the jump size increasing 
with number of monomers N . This value of A corre- 
sponds to the collapse transition of the condensed poly- 
electrolyte (see inset of Fig. [3]). 

We also study the relative fluctuations Xb of the num- 
ber of non-bonded neighbours, where 

N [{np - {ubf] 

It has a peak around A' « 0.89 (see Fig. [5]). The in- 
creasing peak height with the number of monomers iV, 
indicates a divergence in the thermodynamic limit. This 
critical value A' corresponds to the onset of condensation 
of the counterions. At this value of A, the polyelectrolyte 
chain is fully extended and hence this condensation is 
reminiscent of Manning condensation. The data for dif- 
ferent values of N can be collapsed using the scaling form 

Xb « N'^'f [{A - A')N'^-] , iV » 1, (8) 

where 0i and (/)2 are scaling exponents. Data collapse is 
seen for t/fi « 0.20 and 02 ~ 0.15 (see inset of Fig. [5]). 



FIG. 4: The number of non bonded neighbours per monomer 
(ub) as a function of A for different A'^ 
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FIG. 5: The relative fluctuation Xft of the number of non 
bonded neighbours, as defined in Eq. ((7)l, as a function of 
A for different A^. Inset: Data collapse when Xb and A are 
scaled as in Eq. ©, where A' = 1.08 for TV = 50 and A' = 0.89 
otherwise. 

Divergence of Xb with data collapse is a strong indica- 
tion of condensation being a continuous phase transition. 
Earlier discussion of the order of the Manning like con- 
densation has been restricted only to model systems of 
cylinder in three dimensions and disc in two dimensions, 
and a similar continuous transition has been reported in 
such systems [26l - [28j . 

D. Transition from extended phase to collapsed 
phase 

To quantify the transition from the extended phase 
to the collapsed phase, the electrostatic energy per 
monomer Ec and its fluctuations are calculated. Only 
monomer- monomer pairs are used for calculating these 
quantities since we are interested in the configuration of 
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FIG. 6: Mean Coulomb energy per monomer {Ec) as a func- FIG. 7: Fluctuation in the Coulomb energy per monomer Xc 
tion oi A for different N. Inset: {Ec) ~ N'^^^ in the collapsed as a function of A. Inset: Time series of the electrostatic 
phase. energy per monomer, for A'' = 200, and A — 6.50. 



the polyelectrolyte chain. The relative fluctuation Xc in 
the electrostatic energy is defined as 

The mean electrostatic energy per monomer {Ec) in- 
creases with A and has an abrupt jump at Ac ~ 6.25 (see 
Fig. [6|). The jump is more pronounced for higher values 
of N. In the collapsed phase, the shape of the chain is 
roughly spherical and hence (Ec) should scale as N^^-^ 
[38l |. For A > Ac, we confirm that (Ec) scales as N^/^ 
(see inset of Fig. The scaling of (Ec) as N'^/^ along 
with the scaling of Rg as N^^^ (see inset of Fig.[3|) in this 
regime confirm that the polymer is indeed in a collapsed 
configuration. 

The variation of the relative fiuctuation Xc with A is 
shown in Fig.[71 Xc has a peak at Ac, close to the value 
of A at which (Ec) has a discontinuity. The peaks are 
not resolved to the desired accuracy. This is because, 
close to the transition point, the polyelectrolyte chain 
fluctuates between the extended and the collapsed phases 
during the time evolution. The inset of Fig. [7] shows 
a sample time series of Ec, near the transition point, 
for N = 200, where Ec fluctuates roughly between two 
values. This behaviour of Xc coupled with the sharp rise 
in (Ec) (see Fig. [6]) suggests a first order phase transition. 
Our findings are consistent with the previous theoretical 
results 39] that the extended to collapse transition of a 
polyelectrolyte chain is first order. 

Similar divergence in the fluctuation of the electro- 
static energy per monomer has been observed for the 
transition of a polyelectrolyte chain, in a poor solvent, 
from the initial globular phase to the dumbbell and trim- 
bell phases in the absence of explicit counterions jTi] . 



E. Existence of Sausage phase 

The sausage phase was defined in Ref. [1^1 as a col- 
lapsed phase where mean asphericity is non-zero. We 
show below that the transition associated with change 
in asphericity coincides with the collapse transition. We 
define asphericity Y as 

/ X A2+A3 \ 

where A 1.2, 3 are the eigenvalues of the moment of inertia 
tensor with Ai being the largest eigenvalue. The moment 
of inertia tensor G is 

1 ^ 

Ga^J = ^^HaH^, (11) 

4=1 

where Via is the a*'' component of the position vector 
fi. Asphericity Y is zero for a sphere (collapsed globule) 
and one for a linear rod (extended configuration). For 
all other configurations, it has a value between zero and 
one. 

The variation of asphericity with A for different N is 
shown in Fig. [5] For very small values of A, aspheric- 
ity increases corresponding to the extension of the initial 
collapsed phase. In the extended region 0.89 ^ A< 6.25, 
asphericity increases with N and tends to one for large 
N. For A > 6.25, asphericity decreases to zero with N as 
a power law (see inset of Fig. [8]). Thus, for large N, as- 
phericity jumps from one to zero as A crosses Ac ~ 6.25. 
The value of the transition point coincides with the ex- 
tended to collapsed transition discussed in Sec. IIII Dl As 
there is no second transition in asphericity or in (Ec), we 
conclude that the sausage phase suggested in Ref. [15] is 
identical with the collapsed phase and is not a different 
phase. 
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A 

FIG. 8: Asphericity Y (as defined in Eq. (|10|l as a function of 
A. Inset: Variation of asphericity Y with chain length A'^ for 
A — 10.93. The slope of the straight line is —0.33. 

IV. CONCLUSIONS 

Using molecular dynamics simulations, we studied the 
phase diagram and phase transitions between the phases 
of a single polyelectrolyte chain in a poor solvent. The 
counterions were taken care of explicitly while the solvent 
was implicit. The dependence of various physical quan- 
tities on A, a dimensionless number that parametrised 
the strength of the electrostatic interaction, were calcu- 
lated for different chain lengths N. Our main results are 
summarised below. 

We quantified the transition associated with the con- 
densation of counterions on the extended chain. The 
fluctuations of the non-bonded neighbours Xb diverges 



at A' w 0.89 with exponents that indicate a continuous 
transition. While the value of the exponents we obtained 
are not accurate, they indicate that data for different 
system sizes can be explained by one scaling function. It 
would be interesting to study the number of non-bonded 
neighbours and its fluctuations for the good solvent prob- 
lem as well as in the presence of salt. 

We also analysed the sausage phase introduced in 
Ref. [iHl as a new intermediate phase between the ex- 
tended and collapsed phases. We show that the transi- 
tion associated with the discontinuity in asphericity and 
that associated with the extended to collapsed transi- 
tion occur, within numerical error, at the same value of 
A {Ac w 6.25). This, in conjunction with the fact that 
there is no second transition for either quantity, strongly 
suggests that the sausage phase does not exist indepen- 
dent of the collapsed phase, and was probably an artefact 
of earlier simulations [Tsj of chains of small size. 

The extended to collapsed transition of the condensed 
polymer chain was also quantified. The abrupt jump 
m Rg/N^/^, in the mean electrostatic energy {Ec), and 
mean number of non-bonded nearest neighbours (nf,), 
transition of the chain between extended phase and col- 
lapsed phase near the transition point, combined with 
the peaks in the fluctuations of the electrostatic energy 
Xc are evidences for a first order transition. This is con- 
sistent with the theoretical prediction in Ref. [39j . 

In addition, we find that the fraction of condensed 
counterions increases with N for fixed A in the extended 
phase. By studying its dependence on tram and the de- 
pendence of Rg/N on A^, we argued that the increase in 
the fraction is related to longer chains having a smaller 
relative extension. 
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